A neural network interatomic potential for the phase change material GeTe 
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GeTe is a prototypical phase change material of high interest for applications in optical and elec- 
tronic non-volatile memories. We present an interatomic potential for the bulk phases of GeTe, 
which is created using a neural network (NN) representation of the potential-energy surface ob- 
tained from reference calculations based on density functional theory. It is demonstrated that the 
NN potential provides a close to ab initio quality description of a number of properties of liquid, 
crystalline and amorphous GeTe. The availability of a reliable classical potential allows addressing 
a number of issues of interest for the technological applications of phase change materials, which 
are presently beyond the capability of first principles molecular dynamics simulations. 
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I. INTRODUCTION 

Phase-change materials based on chalcogenide alloys 
are attracting an increasing interest worldwide due to 
their ability to undergo reversible and fast transitions be- 
tween the amorphous and crystalline phases upon heat- 
ingi. This property is exploited in rewriteable opti- 
cal media (CD, DVD, Blu-Ray Discs) and electronic 
nonvolatile memories (NVM), which are based on the 
strong optical and electronic contrast between the two 
phases 2 .* 3 -. The material of choice for NVM applications 
is the ternary compound Ge2Sb2Tes (GST). However, 
the related binary alloy GeTe has also been thoroughly 
investigated because of its higher crystallization temper- 
ature and better data retention at high temperature with 
respect to GST. 

In the last few years atomistic simulations based on 
density functional theory (DFT) have provided useful in- 
sights into the properties of phase change materials^—. 
However, several key issues such as the thermal con- 
ductivity at the nanoscale, the crystallization dynamics, 
and the properties of the crystalline/amorphous inter- 
face, just to name a few, are presently beyond the reach 
of ab initio simulations due to the high computational 
costs. The development of reliable classical interatomic 
potentials is a possible route to overcome the limitations 
in system size and time scale of ab initio molecular dy- 
namics. However, traditional approaches based on the 
fitting of comparably simple functional forms for the in- 
teratomic potentials are very challenging due to the com- 
plexity of the chemical bonding in the crystalline and 
amorphous phases revealed by the ab initio simulations. 
A possible solution has been proposed recently by Behler 
and Parrinello 9 ., who developed high-dimensional inter- 
atomic potentials with close to ab initio accuracy employ- 
ing artificial neural networks (NN) . To date, potentials of 
this type have been reported for silico n 10 ! 11 , carbo n 12 ' 13 , 
sodium^, zinc oxidei^ and copper— by fitting large ab 
initio databases. 

Here we describe the development of a classical inter- 
atomic potential for the bulk phases of GeTe employing 



this NN technique. The potential is validated by com- 
paring results on the structural and dynamical properties 
of liquid, amorphous and crystalline GeTe derived from 
NN-based simulations with the ab initio data obtained in 
our previous workil. 



II. METHODS 
A. The Neural Network Method 

Artificial neural networks constitute a class of algo- 
rithms inspired by the properties of biological neural net- 
works, which are widely applied in many different fields 
ranging from weather forecasting to robotics^. In the 
last few years, NN have also been used as a tool to con- 
struct interatomic potential o 19 ' 20 To this aim, the NN are 
exploited as a non-linear technique that allows fitting any 
real-valued function to arbitrary accuracy, without any 
previous knowledge about the functional form of the un- 
derlying proble m 21 ' 22 . In our case, this function is the 
potential energy surface (PES) of the atomistic system 
and the goal of the NN is to construct a functional rela- 
tion between the energy and the atomic configuration. 

Given a number of atomic configurations, for which the 
ab initio energies are known, the NN parameters are ad- 
justed to reproduce these energies in the training process. 
Overfitting, i.e., obtaining a good fit of the training data, 
but performing less accurately when making predictions 
is controlled by testing the performance of the NN for 
an independent test set not used in the fitting. Once 
trained, the NN performs an interpolation to construct 
the potential energy for new atomic configurations at a 
low computational load, which allows performing long 
molecular dynamics runs for large systems. 

The NN methodology overcomes many problems as- 
sociated with classical potentials. First and foremost, 
NN completely obviate the problem of guessing a compli- 
cated functional form of the interatomic potential. This 
form is determined automatically by the NN. Moreover, 
the entire training procedure is fully automated so that 
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FIG. 1. Normalized cumulative histograms of the absolute NN errors in training and test data sets for the energies (a) and 
forces (b). Standard histograms for the same data are presented as insets. 



NN can be readily extended to new regions of the PES. 
Finally, the accurate mapping of ab initio energies en- 
sures that all properties determined by the topology of 
the PES are described with an accuracy comparable with 
that of first principles calculations. NN have been suc- 
cessfully used to interpolate the PES of simple chemical 
systems like small molecules^— and molecule-surface 
interactions^— . 

Until recently, the main limitation of most NN ap- 
proaches has been the small number of degrees of freedom 
that could be described, which confined their application 
to very small systems. This limitation has been over- 
come by Behler and Parrinello who introduced a NN po- 
tential, in which the total energy is expressed as a sum 
of atomic energy contributions depending on the local 
environments 9 . In this approach the energy contribution 
of each atom is evaluated using an individual NN instead 
of using one NN for the total energy of the system. The 
local environment of a given atom is described by a set 
of local parameters called symmetry functions, which in- 
clude radial and angular many-body terms and depend 
on the positions of all neighbors within a specified cut- 
off radius^. The use of symmetry functions instead of 
Cartesian coordinates as NN inputs and the partitioning 
of the total energy into atomic contributions ensures that 
all quantities computed with the NN potential, such as 
energies, analytical forces, and stress tensor are invari- 
ant to translations, rotations, and atoms exchange. Fur- 
thermore, once the fit is obtained, the NN potential can 
be applied to systems containing an arbitrary number of 
atoms. The validity of this approach has been demon- 
strated by reproducing the high-pressure phase diagram 
of silicon, carbon, and sodiu m 10 i 12 i 14 . This scheme has 
been recently extended to binary systems by including 
long range Coulomb interaction between environment de- 
pendent ionic charge o 15 i 35 . 

In the case of GeTe, we are faced with the problem 
of developing a potential suitable to describe both the 
semiconducting crystalline and amorphous phases as well 
as the metallic liquid. As a first step toward the devel- 
opment of a NN potential for GeTe, we neglect long- 



range Coulomb interactions for atoms being separated 
by a larger distance than the cutoff radius of the symme- 
try functions. The resulting "short-ranged" NN poten- 
tial just consists of atomic energy contributions arising 
from the local chemical environments, but it is important 
to note that also short-ranged electrostatic interactions 
are fully taken into account implicitly. More details of 
the employed method can be found elsewhere 9 -^.. The 
resulting potential will not be able to describe the dielec- 
tric response and LO-TO splitting of crystalline GeTe in 
its ferroelectric phase, but it will be suitable to repro- 
duce structural properties of the liquid, amorphous and 
crystalline phases and the dynamical properties of the 
disordered phases. 



B. The NN potential for GeTe 

To generate the NN potential, we fitted the total en- 
ergies of about 30.000 configurations of 64-, 96- and 216- 
atom supercells computed by DFT. We started with a 
relatively small dataset of about 5.000 structures, where 
we considered crystalline configurations, snapshots of the 
liquid phase and of the amorphous phase generated by 
quenching from the melt at ambient conditions and at 
different pressures up to 50 GPa. We also considered 
mixed crystalline/amorphous models generated by par- 
tially crystallizing the amorphous phase by means of the 
metadynamics technique- - . All these configurations were 
generated within ab initio molecular dynamics simula- 
tions at different temperatures (up to 3000 K) with the 
code CP2K— . This first dataset was then expanded 
by adding randomly distorted structures of the initial 
dataset at slightly different pressures and temperatures 
and models with slight deviations from the perfect stoi- 
chiometry. The refinement of the potential was achieved 
by inserting in the data set the ab initio energy of config- 
urations generated by Molecular Dynamics simulations 
(see below) using the not yet refined NN potential. 

Total energy convergence has to be strictly guaran- 
teed for every point of the dataset to make the NN input 
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fully consistent. Therefore, we had to perform Brillouin 
Zone (BZ) integration over a dense 4x4x4 Monkhorst- 
Pack (MP)2& mesh for the 64- atom cell and employing 
meshes of a corresponding k-point density for the larger 
systems. The ab initio energies were calculated using the 
QUANTUM-ESPRESSO package^. The Perdew-Burke- 
Ernzerhof (PBE)i£ exchange-correlation functional and 
norm conserving pseudopotentials were employed, con- 
sidering only the outermost s and p electrons in the va- 
lence shell. The Kohn-Sham orbitals were expanded in a 
plane waves basis up to a kinetic energy cutoff of 40 Ry. 
These settings ensure convergence of the total energy to 
2 meV/atom. 

The best NN fit we found employs three hidden lay- 
ers with 20 nodes each. Sigmoidal activation functions 
were used in the nodes of the hidden layers, while a linear 
function was used for the output node. Details concern- 
ing the architecture of the NN can be found elsewhere^. 
The local environment of each atom is defined by the 
value of 159 symmetry functions (see Ref.— for details) 
defined in terms of the positions of all neighbors within a 
distance cutoff of 6.88 A. The generation of the NN po- 
tential and the calculation of the forces for the MD simu- 
lations were performed with the NN code RuNNer 41 . We 
used the TINKER^ code as MD driver. The time step 
for the MD runs was set to 0.2 fs, and constant temper- 
ature simulations were performed using the Berendscn 
thermostat^. 

The results of the fitting process of the NN potential 
are summarized in Fig. [1] The root mean square error 
(RMSE) for the energy is 5.01 and 5.60 meV/atom for 
the training and the test set, respectively, while the RM- 
SEs of the forces are 0.46 and 0.47 eV/A, for the two 
sets. Among all structures considered, only a negligi- 
ble fraction shows noticeable absolute errors, up to 25.8 
meV/atom and 11.2 eV/A for energies and forces (see 
Fig. [T] insets). We have found that these configurations 
correspond to high-energy structures that are not visited 
in MD simulations carried out in the present work. 



III. RESULTS 

A. Crystalline phase 

The equilibrium geometry of the trigonal phase of crys- 
talline GeTe (i?3m space group) was obtained by opti- 
mizing all structural parameters consisting of the lattice 
parameter a, the trigonal angle a and the internal pa- 
rameter x that assigns the positions of the two atoms 
in the unit cell, namely Ge at (x,x,x) and Te at (-af,- 
x,-x)^. The residual anisotropy in the stress tensor at 
the optimized lattice parameter at each volume is below 
0.02 kbar. The energy versus volume data were fitted 
with a Murnaghan equation of state 4 ^. The theoretical 
structural parameters of the trigonal phase of GeTe at 
equilibrium are compared in Tab. U with experimental 
data 4 ^ and DFT results obtained with a 12 x 12 x 12 



MP k-point mesh in the BZ integration. DFT data are 
similar to those reported previously^. The length of the 
short and long Ge-Te bonds are also given. The struc- 
ture of trigonal GeTe can be seen as a distorted rocksalt 
geometry with an elongation of the cube diagonal along 
the [111] direction and an off-center displacement of the 
inner Te atom along the [111] direction, which moves to 
a distance d from the Ge atom at the vertex as shown 
in Fig. [2^,. The energy gained by the off-center displace- 
ment is analyzed by varying the distance d at fixed lattice 
parameters a and a . The resulting energy as a function 
of d is reported in Fig. [5b for the NN and the DFT calcu- 
lations. We note that the DFT values were not included 
in the training set but were recalculated for investigating 
the quality of the NN potential only. The double well 
potential identifies the two possible ferroelectric configu- 
rations while the maximum corresponds to a paraelectric 
configuration. 



Structural parameters 


NN 


DFT 


Exp. 


a (A) 


4.47 


4.33 


4.31 


a 


55.07° 


58.14° 


57.9° 


Volume (A 3 ) 


55.95 


54.98 


53.88 


X 


0.2324 


0.2358 


0.2366 


Short, long bonds (A) 


2.81, 3.31 


2.85, 3.21 


2.84, 3.17 



TABLE I. Structural parameters of the trigonal phase of crys- 
talline GeTe from NN and DFT calculations and from the ex- 
perimental data of Ref. 0. The lengths of the short and long 
bonds are also given. 
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FIG. 2. a) Trigonal structure of crystalline GeTe, which can 
be seen as a cubic rocksalt geometry with an elongation of 
the [111] diagonal. Ge atoms are on the vertexes of the cell 
while Te atom stands in the center. The distance between Ge 
and Te atoms along the [111] direction is indicated by d. b) 
Energy of trigonal GeTe as a function of d at fixed values of 
the lattice parameters at the theoretical equilibrium geometry 
(cf. Table 15. 
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As a further validation of the potential, we computed 
the difference in energy between the trigonal phase and 
an ideal rocksalt phase at their equilibrium volumes at 
zero temperature that amounts to 44 meV/atom or 55 
meV/atom in NN and DFT calculations, respectively. 





Cll 


Cl2 


C13 


Cl4 


C33 


C44 


DFT 


92 


18 


22 


35 


40 


24 


NN 


73 


10 


30 


24 


36 


20 



TABLE II. Elastic constants (GPa) of trigonal GeTe from 
DFT and NN calculations. 




The elastic properties of trigonal GeTe were inves- 
tigated by computing the elastic constants from finite 
deformations of the lattice parameters. The NN and 
DFT results are compared in Table [TTJ The elastic con- 
stants obtained here with the PBE functional are some- 
how softer than those obtained with the LDA functional 
in Ref. H3- The bulk modulus obtained either from the 
elastic constants or from the equation of state is 34 GPa 
and 33 GPa for the NN and DFT calculations. 



B. Liquid phase 

The liquid phase of GeTe was simulated by a 4096- 
atom model at 1150 K. Total and partial pair correlation 
functions are compared in Fig. [3] with results from our 
previous ab initio simulations in a small 216-atom cell at 
the same temperature^. Results from the NN simula- 
tions of a 216-atom cell are also reported. The density 
of 0.03156 atoms/A 3 is the same for all simulations and 
corresponds to the value chosen in the ab initio simula- 
tions of Ref. Qj], which is close to the experimental density 
of the amorphous phase 48 . Distributions of coordination 
numbers are reported in Fig. 2] as computed by integrat- 
ing the partial pair correlation functions up to the cut- 
off shown in Fig. [3] Average coordination numbers are 
given in Table IIIII while angle distribution functions are 
shown in Fig. [SJ The agreement between NN and ab ini- 
tio data is excellent. The NN results obtained with 216- 
atom and 4096-atom cells are extremely similar, which 
demonstrates that structural properties of the liquid can 
be reliably described by the cells few hundred atoms large 
we used in our previous ab initio works^ 7 - ' 17 ' 49 ' 50 . The 
self-diffusion coefficients computed from NN simulations 
are also in good agreement with the ab initio results of 
Ref. [H as shown in Table IIVI These latter data refer to 
simulations at 1000 K to enable a comparison with pre- 
vious ab initio results obtained at this temperature. 



FIG. 3. Total and partial pair correlation functions of liquid 
GeTe from a NN molecular dynamics simulation at 1150 K 
with a 4096-atom and a 216-atom cell, compared with results 
from an ab initio simulation at the same temperature using 
a 216-atom cell^. NN results are obtained by averaging over 
a NVE run 40 ps long at the average temperature of 1150 K. 
The vertical lines are the interatomic distance threshold used 
to define the coordination numbers, 3.0 A, 3.22 A, and 3.0 
A for Ge-Ge, Ge-Te and Te-Te bonds, respectively. 
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FIG. 4. Distributions of coordination numbers of Ge and 
Te atoms in liquid GeTe at 1150 K. Results from NN (4096- 
atom and 216-atom) and ab initio (216-atom)^ 7 - simulations 
are compared. 



C. Amorphous phase 

The structural properties of a-GcTe and a-GST have 
been elucidated recently by ab initio simulations^ " 8 ' 17 ' 51 . 
In these systems Ge and Te atoms are mostly four- 
coordinated and three-coordinated, respectively. Te 
atoms are in a defective octahedral-like environment, 
which resembles the local environment of the correspond- 
ing crystalline phases. The majority of Ge atoms are in 
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With Ge 


With Te 


Total 




DFT 


NN216 


NN40U6 


DFT 


NNaie 


NN 4 096 


DFT 




NN 40 96 


Ge 


1.00 


1.11 


1.15 


2.71 


2.78 


2.67 


3.71 


3.89 


3.82 


Te 


2.71 


2.78 


2.67 


0.26 


0.28 


0.26 


2.97 


3.07 


2.93 



TABLE III. Average coordination numbers for different pairs 
of atoms computed from the partial pair correlation functions 
of liquid GeTe from a NN molecular dynamics simulation at 
1150 K with a 4096-atom and a 216-atom cell (cf. Fig. [3|, 
compared with results from an ab initio (DFT) simulations of 
a 216-atom cell at the same temperature—. The interatomic 
distance thresholds defined in Fig. [3] are used. 




20 40 60 80 100 120 140 160 180 
(deg) 

FIG. 5. Total and partial angle distribution functions of liquid 
GeTe from a NN molecular dynamics simulation at 1150 K 
with a 4096-atom and a 216-atom cell, compared with results 
from an ab initio simulation at the same temperature using 
a 216-atom cell^. Partial distributions refer to X-Ge-Y and 
X-Te-Y triplets (X,Y= Ge or Te). 





NN 


DFT 


D Ge (10~ 5 cm 2 /s) 


4.96 


4.65 


D Te (10 -5 cm 2 /s) 


3.62 


3.93 



TABLE IV. Diffusion coefficient of Ge and Te atoms in the 
4096-atom model of liquid GeTe at 1000 K. NN results were 
obtained from the slope of the mean square displacement ver- 
sus time. The same values within two decimals are obtained 
from the integral of the velocity- velocity autocorrelation func- 
tion. DFT data of a 216-atom cell at the same temperature 
are from Ref. 0. 



a defective octahedral environment too, but about one 
quarter of Ge atoms are in a tetrahedral-like geometry. 
The presence of homopolar Ge-Ge (and, in the case of 
GST, Ge-Sb) bonds favors the tetrahedral coordination. 

In the following, we compare the structural properties 
of models of amorphous GeTe generated by NN and ab 
initio simulations. The NN amorphous phase was gen- 
erated by quenching the molten sample from 1150 K to 
room temperature in 100 ps. Average properties are ob- 
tained from a NVE simulation 40 ps long at an aver- 
age temperature of 300 K. Doubling or even tripling the 
quenching time (up to 300 ps) does not introduce siz- 
able changes in the structural and vibrational properties 
of our NN model of amorphous GeTe. Structural prop- 
erties are described in Figs. [BJIHl The partial pair cor- 
relation functions of our NN models are compared with 
ab initio data in Fig. [6] By decreasing the system size 
from 4096-atom to 1728- atom one obtains essentially the 
same results. By using instead a small 216-atom cell fluc- 
tuations in properties related to Ge-Ge homopolar bonds 
are found by looking at ten different independent mod- 
els. The quantities averaged over ten models are reported 
hereafter and are very close to those of the larger models. 
The fluctuations among the different 216-atom models 
can be appreciated in Figs. S1-S2-S3 of the additional 
material (EPAPS^). 

The distribution of coordination numbers and their av- 
erage values are reported in Fig. [7] and Table [V] for NN 
and ab initio^ simulations. Bond angles distribution 
functions are reported in Fig. [5] 
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FIG. 6. Total and partial pair correlation functions of amor- 
phous GeTe from a NN molecular dynamics simulation at 
300 K with a 4096-atom and a 216-atom cell, compared with 
results from an ab initio simulation at the same temperature 
using a 216-atom cell—. The vertical lines are the interatomic 
distance thresholds used to define the coordination numbers, 
3.0 A 3.22 A and 3.0 A for Ge-Ge, Ge-Te and Te-Te bonds, 
respectively. 

The agreement between NN and DFT data is over- 
all very good. The largest discrepancy is on the height 
of the first peak of the Ge-Ge pair correlation function. 
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FIG. 7. Distribution of coordination numbers of Ge and Te 
atoms in a-GeTe from NN and DFT— simulations with 4096- 
and 216-atom cells. The interatomic distance thresholds de- 
fined in Fig. [6] are used. 






With Ge 


With Te 


Total 




DFT 


NN 2 16 


NN4096 


DFT 


NN 21fj 


NN4096 


DFT 


NN 21B 


NN4096 


Ge 


0.76 


0.78 


0.88 


3.24 


3.31 


3.22 


4.00 


4.09 


4.10 


Te 


3.24 


3.31 


3.22 


0.02 


0.04 


0.05 


3.27 


3.35 


3.27 



TABLE V. Average coordination number for different pairs of 
atoms computed from the partial pair correlation functions of 
amorphous GeTe from NN molecular dynamics simulations at 
300 K with a 4096-atom and a 2f6-atom cell, compared with 
results from DFT simulations at the same temperature and a 
216-atom celli£. The interatomic distance thresholds defined 
in Fig. [6] are used. 



However, this misfit is partially due to fluctuations in 
the number of Ge-Ge bonds in the still small 216-atom 
cell (cf. Figs. S1-S2 in EPAPS^). In fact, by generat- 
ing ten different 216-atom models with the NN potential 
we observed sizable fluctuations in the average Ge-Ge 
coordination number ranging from 0.6 to 0.9. A more 
compelling comparison would require the availability of 
several independent ab initio models of a-GeTe. 

Another discrepancy with the DFT results is the pres- 
ence of a small peak at around 60° in the angle distribu- 
tion function due to a very small fraction of three mem- 
bered rings (see below). Following our previous works^, 
we quantified the fraction of Ge atoms in a tetrahe- 
dral geometry by computing the local order parameter 
q = 1 — | 'Yli > k{\-\-cos0ijk)' 2 where the sum runs over the 
pairs of atoms bonded to a central atom j. q = 1 for the 
ideal tetrahedral geometry, q = for the six-coordinated 
octahedral site, and q — 5/8 for a four-coordinated de- 
fective octahedral site. The distribution of the local or- 
der parameter q for Ge atoms is reported in Fig. [9] for 
different coordination numbers. The q distribution for 4- 
coordinated Ge is bimodal with peaks corresponding to 
defective octahedra and tetrahedra. In contrast, the q- 



I 1 lJ-=C ' ' ' ' =^ 

20 40 60 80 100 120 140 160 180 

8 (deg) 

FIG. 8. Total and partial angle distribution functions of amor- 
phous GeTe from a NN molecular dynamics simulation at 300 
K with a 4096-atom and a 216-atom cell, compared with re- 
sults from an ab initio simulation at the same temperature 
and a 216-atom celtf. Partial distributions refer to X-Ge-Y 
and X-Te-Y triplets (X,Y= Ge or Te). 



distribution for Te does not show any signature of the 
tetrahedral geometry (cf. Fig. [9]). We estimated the 
fraction of tetrahedral Ge atoms by integrating the q- 
distribution of 4-coordinated Ge from 0.8 to 1. This pro- 
cedure was demonstrated to provide reliable values for 
the fraction of tetrahedral Ge from the analysis of the 
Wannicr functions that allow a direct identification of 
the tetrahedral geometry in terms of the electronic struc- 
tur o 50 ' 51 . In fact, Ge in tetrahedral sites has four bonding 
sp 3 -like Wannier functions, while Ge in defective octahe- 
dra has three p-like bonding Wannier functions and one 
s-likc lone pair. The fraction of tetrahedral Ge atoms for 
the ten 216-atom NN models, for the large 4096-atom 
model and for the ab initio 216-atom model of Ref. H3 is 
given in Fig. 1101 By chance the ab initio value is very 
close to the value obtained by averaging over the ten 216- 
atom NN models while the fraction of tetrahedral Ge in 
the large 4096-atom cell is 24 %, a value very close to the 
average over the ten small cells and also to the result ob- 
tained for a 1728-atom model (22 %). The concentration 
of tetrahedra depends on the fraction of Ge-Ge homopo- 
lar bonds, which is the property more affected by finite 
size effects in the small 216-atom cell. 

Turning now to the medium range order, we report in 
Fig. [H] the distribution of ring lengths computed accord- 
ing to Ref. [53| for the large and small NN models and for 
the small DFT model of Ref. Il7j. The amorphous phases 
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FIG. 9. Order parameter q for tetrahedricity (see text) for 3-, 4-, and 5-fold coordinated Ge and Te atoms in amorphous GeTe 
from NN molecular dynamics simulations at 300 K with a 4096-atom and a 216-atom cell, compared with results from an ab 
initio simulation at the same temperature employing a 216-atom cell^. 




3 4 5 6 7 8 9 10 4 6 8 10 12 14 1 6 1 8 20 

Model number Ring length (n. of atoms) 



FIG. 10. Fraction of tetrahedral Ge atoms in amorphous GeTe 
from ten statistically independent 216-atom models generated 
with the NN potential, compared with results from a single 
ab initio 216-atom model^. The average (AVG) over the ten 
NN models and from the 4096-atom model are also reported. 



of GeTe and GST have been shown to display a large 
concentration of nanocavities^. The distribution of the 
volume of nanocavities computed according to the defi- 
nition of Ref. [Hi] and the algorithm of Ref. [55ll56l is com- 
pared in Fig. [12] for the NN and DFT models. The same 
scheme for the calculation of nanocavities was applied in 
our previous works on different phase change materials^ 9 . 
These comparisons show that the agreement between NN 
and DFT results is very good for the medium range order 
as well. 

We also optimized the density of the amorphous model 
at zero temperature by interpolating the energy-volume 
points with a Murnaghan equation of state. We obtained 
an equilibrium density of 0.03351 atoms A~ 3 to be com- 
pared with the value of 0.03156 atoms A -3 resulting 
from the ab initio equation of state of a 216-atom cell 
with the BZ integration restricted to the r-point^. The 
experimental equilibrium density^ of a-GeTe is 0.03327 
atoms A~ 3 . The NN and ab initio bulk moduli of a-GeTe 



FIG. 11. Distribution of ring lengths in the 4096- and 216- 
atom NN models and in the DFT 216-atom model of Ref.Qj]. 




2 4 6 8 10 12 14 16 18 20 
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FIG. 12. Distribution of the volume of cavities in the 4096- 
and 216-atom NN models and in the DFT 216-atom model of 
Refs. [HM 



are 17 GPa and 14 GPa, respectively. 

Concerning the vibrational properties, the phonon den- 
sity of states of amorphous GeTe from NN and ab initio 
simulations are compared in Fig. 1131 Phonon frequencies 
are computed by diagonalizing the dynamical matrix ob- 
tained in turn from the variation of atomic forces due to 
finite atomic displacements 0.02 A large. Only phonons 
with the periodicity of our supercells (r-point phonons) 
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were considered. Ab initio phonons are computed in a 
216-atom cellii while NN phonons are obtained from the 
4096-atom and 216-atom cells. Projections on the differ- 
ent type of atoms (Te, Ge in tetrahedral and defective 
octahedral geometries) are also shown. 

In an amorphous material, phonons display localiza- 
tion properties, which depend on frequency. To address 
this issue and following our previous ab initio workai£, we 
computed the inverse participation ratio (I PR) of the j- 
th vibrational mode defined as 



IPR 



|eQ»|- 



(1) 



FIG. 13. Phonon density of states of amorphous GeTe from 
4096- and 216-atom NN models. Results for a single 216- 
atom DFT mode^ are also reported. Only phonons at the 
supercell T point are considered. The higher smoothness of 
4096-atom NN data is due to the larger supercell used. Pro- 
jections of the DOS on the different atomic species (Te atoms 
and Ge atoms in tetrahedral and defective octahedral geome- 
tries) are also shown. 




50 100 150 200 250 300 

to (cm -1 ) 



Here e(j, k) are phonon eigenvectors and the sum over 
k runs over the N atoms in the unit cell with masses 
M K . According to this definition, the value of IPR varies 
from 1/N for a completely delocalized phonon, to one for 
a mode completely localized on a single atom. The val- 
ues of IPR for the NN and ab initio models of a-GeTe 
are reported in Fig. [TJ] The NN potential reproduces 
the strong localization on tetrahedra of phonons above 
200 cm -1 . The overall shape and frequency range of the 
phonon DOS is reasonably reproduced by the NN poten- 
tial. A discrepancy is present in the relative height of the 
two main structures at 50 cm -1 and 150 cm -1 , which, 
however, might be partially due to the still small size of 
the 216-atom cell. Actually the DOS of the different 216- 
atom models generated by the NN potentials are all very 
similar, but somehow different from the large 4096-atom 
cell due to a size effect. 
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FIG. 14. Inverse participation ratio of phonons in a-GeTe 
from 4096- and 216-atom NN models. Results for a single 
216-atom DFT model^ are also reported. 

IV. CONCLUSIONS 

In summary, a NN potential for the phase change ma- 
terial GeTe has been created and tested to reproduce the 
properties of crystalline, liquid and amorphous phases. 
The NN potential has been validated by comparing the 
results on structural and dynamical properties of the bulk 
phases of GeTe with our previous data from DFT calcu- 
lations 1 !. The development of a classical potential with 
close to ab initio accuracy represents a breakthrough in 
the simulation of phase change materials, as it will allow 
addressing several key issues on the properties of this 
class of materials that are presently beyond the reach of 
ab initio simulations. The study of thermal conductivity 
in the amorphous phased and the dynamics of homo- 
geneous and heterogeneous crystallization are few exam- 
ples of useful follow-on developments of this work, which 
promises to improve our microscopic understanding of 
the operation of phase change memories. 
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Fig. S 15. Total and partial pair correlation functions of amorphous GeTe from NN molecular dynamics 
simulations at 300 K of ten statistically indipendent 216-atom models. Average values are also reported. 
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Fig. S16. Distribution of coordination numbers of Ge and Te atoms of amorphous GeTe from NN molecular 
dynamics simulations at 300 K of ten statistically indipendent 216-atom models. Average values are also 
reported. 
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Fig. S 17. Total and partial angle distribution functions of amorphous GeTe from NN molecular dynamics 
simulations at 300 K of ten statistically indipendent 216-atom models. Average values are also reported. 
Partial distributions refer to X-Ge-Y and X-Te-Y triplets (X,Y= Ge or Te). 



